CoefRcient of restitution of colliding viscoelastic spheres 
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We perform a dimension analysis for colliding viscoelastic spheres to show that the coefficient of 
normal restitution e depends on the impact velocity g as e = 1 — "fig^^^ +72<7^^^ =p . . ., in accordance 
with recent findings. We develop a simple theory to find explicit expressions for coefficients 71 and 
72. Using these and few next expansion coefRcients for e{g) we construct a Pade-approximation 
for this function which may be used for a wide range of impact velocities where the concept of 
the viscoelastic collision is valid. The obtained expression reproduces quite accurately the existing 
experimental dependence e{g) for ice particles. 
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I. INTRODUCTION 

The change of relative velocity of inclastically colliding 
particles can be characterized by the coefficient of resti- 
tution e. The normal component of the relative velocity 
after a collision g' = • e follows from that before the 
collision g = V12 ■ e via 



g 



£5 



(1) 



where vi, V2 and v(, V2 are respectively the veloci- 
ties before and after the collision, while the unit vector 
e = fi2/ \ri2\ gives the direction of the inter-particle vec- 
tor ri2 = ^1 — ^2 at the instant of the collision. 

From experiments as well as from theory it is well 
known that the coefficient of normal restitution e is not a 
constant but it depends sensitively on the impact velocity 
|^-|rT| . Although most of the results in the field of gran- 
ular gases have been derived neglecting this dependence 
but using a velocity-independent coefficient of restitution 
(e.g. [p2|-p8[) it has been shown that the impact- velocity 
dependence of the coefficient of restitution has serious 
consequences for various problems in granular gas dy- 
namics ll-lll. 

The equation of motion for inclastically colliding 3D- 
spheres has been addressed in p^-p^, where the Hertz 
contact law IGTI 



Pi 



3/2 



2Y 



3(1 - 1/2 



(2) 



for the clastic inter-particle force has been extended 
to account for the viscoelasticity of the material which 
causes the dissipative part of the force 



(3) 



Here, ^ is the compression of the particles during the col- 
lision ^ = Ri + R2 — 1^1 — (i?i, i?2 and fi, ^2 are 



the radii and the positions of the spheres), Y and v are 
respectively the Young modulus and the Poisson ratio of 
the particle material, = R1R2/ (Ri + R2), and the 
dissipative parameter A reads 



3 {3m + 27ji) 



(4) 



The viscous constants rji, 772 relate the dissipative stress 
tensor to the deformation rate tensor p5| , p6|j28| . The 

same functional dependence of F^iss (^S,, Cj has been ob- 
tained in using a different approach. 

The equation of motion for inclastically colliding 
spheres reads, therefore, 



p 



(5) 



with e(o) = o, m = g, 



and with m'^^ = mim2/(mi -f- 7712) {mi, 111,2 are the 
masses of the colliding particles) . To obtain the depen- 
dence of the restitution coefficient on the impact veloc- 
ity for 3D-spheres the equation of motion (@) was solved 
numerically |p^ , p^ [2^ and analytically where the 
velocity-dependent restitution coefficient has been ob- 
tained as a series in powers of 5^^^: 



e = 1 



(6) 



C2 -A 



The first coefficients Ci = 1.15344 and C2 = 0.79826 
were evaluated anal ytic ally and then confirmed by nu- 
merical simulations ||32|| . 

Although in |^2| a general method of derivation of all 
coefficients of the expansion (^ has been proposed, to ob- 
tain these, extensive calculations have to be performed. 
This approach does not provide closed-form expressions 



1 



for the coefScients, but rather gives them in terms of con- 
vergent series which are to be evahiated up to the desired 
precision. 

In the present study we show that a dimension anal- 
ysis ahows to obtain the functional form of the e{g)- 
dependence for the elastic and dissipative forces. Within 
the framework of this analysis we reproduce the depen- 
dence (H) up to numerical values of coefficients Ck- A 
similar approach has been used by Tanaka et al. to 
prove that the constant coefficient of restitution is not 
consistent with physical reality (see also ) . We also 

develop a simple approximative theory, which gives a con- 
tinuum fraction representation for e{g) and a closed-form 
expressions for Ci and C2 with the same numerical values 
as above. Using then coefficients Ci, . . . C4 (with C3 and 
C4 evaluated in Appendix in accordance with the general 
scheme of Ref. [^) we construct a Pade-approximation, 
which reproduces fairly well the experimental data for 
colliding ice-particles ||]. 



In the last Eq. (|0|) we supplemented the pre-coUisional 
initial conditions at r = with the after-coUisional con- 
ditions at T — Tc (t is the dimensionless time and Tc is 
the dimensionless duration of the collision) . These follow 
just from the definition of the restitution coefficient. We 
point out that all dependence on the initial impact ve- 
locity on any quantity of the problem, including e (this 
is just the dimensionless after-collisional velocity) comes 
only through the constant 5, which reads 



S{g) = D2 



a 



2Di 



1 + 7 
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Hence, e{g) = e {6{g)). A similar result for e — > 0, /3 
and a = 3/2 has been obtained in 

If we assume that the restitution coefficient does not 
depend on the impact velocity g then follows 



(12) 



II. THE DIMENSIONAL ANALYSIS 

To perform the general dimensional analysis we adopt 
the following form for the elastic and dissipative forces: 



m 



This general form (at least for small ^ and ^) follows from 
the fact that both elastic and dissipative forces vanish at 
^ = and ^ = 0, respectively. With these notations the 
equation of motion for colliding particles reads 

I + DiC +D2Ci^ = 0, 
with e(0)=0, m^g, 



(7) 



where g has already been introduced. Now we choose as 
the characteristic length ^0 of the problem the maximal 
compression for the elastic case. It may be found from 



the condition that the initial kinetic energy m g /2 
equals the maximal elastic energy rrf^ DiS^^"^"^ / {a 4- 
which yields 



1), 



Co 



2Di 
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(8) 



Choosing then the characteristic time of the problem as 
To = Co /5, we construct new dimensionless variables 

C = C/Co, i^i/g, I -A' (9) 



and recast the equation of motion into dimensionless 
form: 
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i + 5{g)Ci^' 
C(0)-0, i(o) = i 



. ^" = with 



(10) 



For a linear dependence of the dissipative force on the ve- 
locity, i.e. for (3 = 1 (this seems to be the most realistic 
for small ^) one obtains a constant restitution coefficient 
for 

• the linear elastic force, Fc\ ^ i.e. a = 1. The 
condition (|l^) implies 7 = 0, and thus the linear 
dissipative force Fdiss ~ 

• the Hertz law for 3D-spheres (||) a = 3/2, there- 
fore 1 = J and i^diss ~ C C^^^ provides a constant 
restitution coefficient. 



We now ask the question: What kind of 6(17) depen- 
dence corresponds to the forces which act during colli- 
sions of viscoelastic particles? It may be generally shown 
p5| , p6 37 1 that the relation 



(13) 



between the dissipative and elastic forces with the dissi- 
pative constant A given in eq. (^) holds, provided three 
conditions are met |Q: 

(i) The elastic part of the stress tensor depends lin- 
early on the deformation tensor p8| |. 

(ii) The dissipative part of the stress tensor depends 
linearly on the deformation rate tensor pS] . 

(iii) The conditions of quasistatic motion are provided, 
i.e. (? <C c, Tvis ^ Tc [ p5|p^ (here c is the speed of 
sound in the material of particles, Tvis is relaxation 
time of viscous processes in its bulk). 
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From this follows that P = 1, -f = a — 1, and thus 
the constant restitution coefficient may be observed only 
for collisions of cubic particles with surfaces normal to 
the direction of collision. We wish to emphasize that this 
conclusion comes from the general analysis of viscoelastic 
collisions. 

Let us discuss now collisions between spheres with elas- 
tic and dissipative forces as given by (^) and (^), respec- 
tively. For these we have rn°^Z?i = p, a = 3/2 and 
TO°*Z?2 = I^P; 7 = 1/2, /3 = 1 which yields the func- 
tional dependence for S{g) and e{g) respectively: 



2 V4 



3/5 



A 



2/5 



,1/5 



= e A 



2/5 



,1/5 



(14) 
(15) 



(skipping the prefactor of S{g) in the last equation) in 
accordance with (^ as found previously. 



III. THE RESTITUTION COEFFICIENT FOR 
SPHERES 



Using 2t = tj^ it is convenient to write the equation 

dEii) 



of motion for a collision in the form 



if5/2 

2^ 



1/2 



m = 0; m = 1 

where we introduce the mechanical energy 



1 

2^ 



if5/2 

2^ 



(16) 



(17) 



To find the energy losses in the first stage of the colli- 
sion, which starts with zero compression and ends in the 
turning point with maximal compression 



^° dE 



-d^ 



Co 



(18) 



one needs to know the dependence of the compression 

rate ^ as a function of the compression ^. 

For the case of elastic collisions, the maximal compres- 
sion is ^0 = Ij according to the definition of our dimen- 

sionless variables. Hence, the dependence follows 
from the conservation of energy: 



1 - e^/' 



(19) 



e(0«Vl-(C/eo)^/' 



(20) 



which also gives vanishing velocity ^ at the turning point 
^0- Integration in Eq. (|l^ ) may be performed yielding 



1 



2^0 



5/2 



1 



- 2 = --^^^0 



3/2 



(21) 



where we take into account that i?(Co) = i E{0) 



and introduce a constant 



i/Vl-a:5/2 =. 



5r(i) 



(22) 



Consider now the inverse collision, which is defined as a 
collision which starts with velocity e g and ends with ve- 
locity g. According to the concept of the inverse collision 
introduced in [ |3^ (which is a useful auxiliary model), 
it is characterized by a negative damping (the energy 
"is pumped" into the system during the collision). The 
maximal compression is the same in both collisions, 
the direct and the inverse. 

Rescaling equation of motion for the inverse collision 
in the very same way as for the direct collision yields 



dm 

di 

m 



1/2 



0, m 



(23) 



This suggests the following approximative relation for 
^(^) during the inverse collision: 



k>^Vi-(^7lo)^/^ 



(24) 



with the additional prefactor e, which is the initial veloc- 
ity in the inverse collision. 

Integration of the energy gain for the first stage of the 
inverse collision (which equals up to its sign the energy 
loss in the second stage of the direct collision H^) may be 
performed just in the same way as for the direct collision, 
yielding the result 



1 



2^° 2 



2^ ■ 



(25) 



Mul- 



where we again use E{^q) — ^S,q and E{0) = 
tiplying Eq. ( pl| ) by e and summing it up with Eq. (|2£^ 
we obtain a simple approximative relation between the 
restitution coefficient and the (dimensionless) maximal 
compression: 



The velocity ^ vanishes at the turning point ^ = 1. For 
inelastic collisions the maximal compression is smaller 
than 1, therefore, one can write an approximation rela- 
tion for the inelastic case: 



(26) 



Substituting this into Eq. ( pT| ) we arrive at an equation 
for the restitution coefficient 
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(27) 



The formal solution to this equation may be written as a 
continuum fraction (which does not diverge in the limit 
g oo): 



l + 2Sd (l + 26d{l + ■ ■ ■) 



2/5 



2/5 



(28) 



Another way of representation of the restitution coef- 
ficient e is a series expansion in terms of S. For practi- 
cal applications it is convenient to return to dimensional 
units. We define the characteristic velocity g* such that 



S = 



1 



2d \g 



9 



1/5 



(29) 



with d beeing defined in Eq. (22). Using the definition 
( p^ ) together with Eq. ( p^ ) we find for the characteristic 
velocity 



(5* 



-1/5 _ r(3/5) /3 

21/552/5 r (21/10) V2 



V^cff 



2/5 



Evaluating the numerical prefactor finally yields 



{9I 



*^"^/^ = 1.15344 ( -A 



2/5 



(30) 



(31) 



Note that the numerical constant 1.15344 has to be equal 
to Ci in Eq. (|). 

With this new notation the restitution coefficient reads 



1 I 9 

9 



, 9 



, 9 



(32) 



with oi — 1, 02 — 3/5 (which are exact values), 03 = 
6/25 = 0.24, a4 = 7/125 = 0.056, • • • (which deviate 



from the correct ones, see below). Comparing ( |3^ ) with 
(^) we conclude that our simple approximative theory re- 
produces exactly the coefficients Ci and C2, which were 
found before using extensive analysis |]3^ . 

We also performed rigorous but elaborated calculations 
according to the general scheme of [p2| to find the exact 
coefficients (details are given in the Ap)pendix) 



= 0.315119, a4 = 0.161167, 



or, respectively, 



Ca = -0.483582 , C4 = 0.285279 . 



(33) 



(34) 



Hence, we observe that while the first two coefficients 
fli = 1 and a2 = 3/5 are correctly obtained from the ap- 
proximative theory, the next approximated coefficients 
as, a4 differ from the exact ones. 

For practical applications such as Molecular Dynam- 
ics simulations, however, the expansion ( ^2|) is of limited 
value, due to its divergence for high impact velocities, 
g — > cxD. According to the velocity distribution function 
there is a certain probability that the relative velocity g 
of colliding particles exceeds the limit of applicability of 
(^). Therefore, we use the obtained coefficients to con- 
struct a Pade-approximation for e{g) which reveals the 
correct limits of the boundary conditions, e(0) = 1 and 
e(oo) = 0. Since the dependence e{g) is expected to be a 
smooth, monotonically decreasing function, we choose a 
"1-4" Pade-approximation 



1 + di 



i+d,{f.y+d,{py+d,{^y+d,{py 

(35) 



Standard analysis 
the coefficients 



fields the coefficients dk in terms of 
3§ (see Table |) 



do ~ a4 — 2a-s — 02 + 3a2 — 1 

1 - a2 + as - 2a4 + (a2 - l)(3a2 - 203)] d^^ = 2.5839 

;(a3 - ci2)(l - 202) - 04] do ^ = 3.5839 

as + al{a2 - 1) - a4(a2 + 1)] dg^ = 2.9839 

a4(a3 - 1) + (aa - a2)(ai - 2a3)] d'^ = 1.1487 

2(03 — a2)(a4 - a2a3) - (04 — a^)^ - 03(03 — a^)] d^^ = 0.3265 



TABLE I. The coefiicients of the Pade formula (|35|) as derived from the coefficients 



Using the characteristic velocity g* — 0.32 cm/s for ice 
as a fitting parameter we could reproduce fairly well the 
experimental dependence of the restitution coefficient of 
ice as a function of the impact velocity g in the whole 
range of g (Fig. |^). The discrepancy with the experi- 
mental data at small g follows from the fact that the ex- 



trapolation expression, e = 0.32/gr'^-^'^^ used in Q has an 
unphysical divergence at g —^ and does not imply the 
fail of the theory for this region. The scattering of the ex- 
perimental data presented in |^ is large for small impact 
velocity according to experimental complications, hence 
the fit formula of H cannot be expected to be accurate 
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enough for too small velocities. Moreover in the region 
of very small velocity possibly other than viscoelastic in- 
teractions influence the collision behavior, e.g. adhesion. 



1 - Bridges et al. (1984) 

2 - Fade approximant 




1.5 2 2.5 3 3.5 
Normal Impact Velocily (cm/s) 

FIG. 1. Dependence of the normal restitution coef- 
ficient on the impact velocity for ice particles. The 
solid line - experimental data of [5], dashed line - the 
Pade-approximation (^) with the constants given in 
the Table and with the characteristic velocity for ice 
g* = 0.32 cm/s. 



IV. CONCLUSION 



V. APPENDIX 



The general method of derivation of the expansion co- 
efficients Ck has been given in |^^. Here, we briefly 
sketch the main lines of derivation, and provide some 
details for the particular case of C3 and C4. Since the 
method of derivation is based on the collection of terms 
with different dependence on the initial velocity g, it is 
convenient to use a scaling, somewhat different from that 
used before for the dimensional analysis. Namely, we 
rescale the time as t' — {p/mcs)^^^g^^^t and the length 
as a; = (p/mcs)^/^^ to recast Eq. (^) into the form pO| 



(36) 



with a = I A(p/moff )2/^ and using all the notations intro- 
duced previously. The initial conditions for the rescaled 
Eq. (|3q ) now read a;(0) = and a;'(0) = g^/^. For sim- 
plicity of notations we will keep in what follows t for the 
rescaled time. As it was shown in the trajectory 
may be expanded in terms of \/t as 



x{t) = foii^/^ ^ jj^^3/2 ^^^2 



(37) 



We developed a dimensional analysis for the inelas- 
tic collision of spherical particles. We could show that 
for 3D-spheres the functional form for e{g) agrees with 
that derived previously |32| using a much more compli- 
cated approach. Using a simple approximative theory we 
found a continuum-fraction representation for e{g) and 
obtained explicit expressions for the coefficients of the se- 
ries expansion of the restitution coefficient in terms of the 
impact velocity. The first two coefficients in this series 
coincide with that found previously by numerical evalu- 
ation. We report also a few next coefficients which we 
have derived within the general approach of a previous 
study Using the first four coefficients of this series 
expansion we constructed a Pade-approximation for e{g). 
It reproduces fairly well the experimental data for collid- 
ing ice particles. The obtained relation for the restitution 
coefficient may be used for a wide range of the impact ve- 
locities, provided that the energy loss during a collision 
is attributed to viscoelasticity and that all the other dis- 
sipative processes (plastic deformation, fragmentation of 
particles) may be ignored. 
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Clearly, both, bi and 63 should be zero to avoid diver- 
gence of velocity and acceleration at t = 0. At the same 
time 62 — g'^^^ and 64 = due to the equation of motion 
at vanishing compression. This yields 



x{t) = g^^H + b^V'/'' + bet^ + brf/^ + , 



(38) 



From ( |3q ) one obtains x'{t) and x"{t) which are to be 
substituted into the equation of motion (|3^). One also 
needs ^/x and a;'^/^; the expansions for these in terms of 
^/i read 



2g2/5 252/i 



2g2/5 



r-f... (39) 



and 



^3/2 ^ ^6/5^3/2 _^ 3^2/5^^^3 ^ ^gV%t'/^ + ... . (40) 



Inserting the expansions for x'{t), x"{t), y/x and cc^/^ 
into (|3^) and collecting the orders of t we obtain 



= ( ^65 + ag'A + Gbet+ (^bj + 1 ] t^^ 



(12 bs + 3ag'/%,) f + (^|69 + lag'/%^ 



.5/2 



(41) 



This suggests the coefficients: 



5 



^^5 

be 



15 ^ 



= 



4 
35 



1 



bs 
69 = 0, 

so that the solution for the trajectory finally reads 



(42) 
(43) 
(44) 

(45) 
(46) 



''max(^max) — 



(50) 



holds at this point. With the above expression for the 
trajectory Eqs. ( ^ , ^9|) , the last Eq. ( |5C| ) is an equation 
to determine tmax, which may be then used to find a;,nax- 
This equation, however is a high-order algebraic equation 
for V^max, which is not generally solvable. On the other 
hand, for the undamped collision ^max equals one-half of 
the collision duration and both quantities of interest 
are known lESl: 



15 



(47) 



In order to get the higher orders, which is conception- 
ally simple but requires extensive calculus, we wrote a 
program |^l[| , which by formula manipulations performs 
exactly the steps we described above and which is able 
to find the trajectory up to any desired order. 

Generally, it is convenient to write the solution as a 
series 

x{t) = (xo{t) + ag^'\^{t) + a^g^">X2{t) -f . . 

(48) 

Here xo(i) is a "zero-order" trajectory, which refers 
to the case of undamped collision, the "first-order" tra- 
jectory, Xi{t)^ accounts for damping in linear (with re- 
spect to a) approximation, the "second-order" trajectory, 
X2{t), corresponds to the next approximation ~ a^, etc. 
Here we give our results for these "n-order" trajectories 
up to n = 3, obtained using the above mentioned pro- 
gram up to the order t^^: 



Xq 



35 175 



22 



17/2 



52 



Xi 



15 



^5/2 ^ 



X2 = TT*' 



15 



38 



70 
937 

' 75075 ' 

,^11/2 + 



104125 



-t 



11 



8017625 



13/2 



238875 
871 



61216 
42639187" 



-X 



808500 



(49) 



43943 



1184627 



-t 



21/2 



2475 13513500 3594591000 

To proceed we need to find the maximal compression 
x-ma-K, which is reached at time imax- The point of max- 
imal compression is a turning point of the trajectory, 
where the velocity is zero. Therefore the condition 



Y2 = T2X'^ 



9"'< ('I 



^1 /// / 



3/5 



2/5 



r(l)r(^: 



(51) 



For a viscoelastic collision imax certainly differs from 
t°/2, so that imax = + ^t. If the dissipation param- 
eter a is not large, the deviation St is presumably small, 

therefore we expand ^'(tinax) — x' {jf ^ ™ terms of 
St: 



g ^ X (^niax) — 



(52) 



-ag 



1/5 



st^'l ( ^ ) + ( ^ 



Stx^ I — 



= 0, 



where we use representation ( |4^ ) for the trajectory. The 
deviation St, vanishes at a = and, thus, suggests the 
expansion in terms of a: 



St = Tia + T2Ci' + T-^O? . 



(53) 



Substituting St, given by Eq. (|53|), into (^2|) and collect- 
ing terms of the same order of a yields 



0, 



(54) 



with the abbreviations 



(55) 
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Vi . r,.'i (I) + (I) + (I) + i'^rW; (I) + (I) + 9'«n4 ) + ^='=4 (| ) 

The conditions Yfe = for /c = 0, . . . 3 together with Eq. ( ^5| ) allows to express the constants ri, T2, T3, etc., in terms 
of functions Xi{t), X2{t), X3{t), etc., and their time derivatives taken at time (i2/2): 




We do not write the expression for T3, since due to the special properties of the problem, i.e. due to the fact that 
^2/2) — the value T3 is not needed for calculation of e up to fourth order of a. The functions xi{t), X2{t), X3{t) 



are known and given by Eqs. (|49|), so that the constants ri and T2 may be found explicitly. 
Writing the maximal compression as 



9 



4/5 



Xq 



ag^^^xi 



ag^/^X2 



X3\^-^+5t 



(57) 



and expanding this in terms of 5t, using then representation of 5t as 5t 
collecting terms of the same order of a we obtain 



C?T2 



•, with Ti, T2 from (pfl) and 



^4/5 



(58) 



where 2/01 ■ • ■ 2/3 are pure numbers: 



yo = ( ^ ) = 1.093362 
= xi I ^ I = -0.504455 



J/2 



2^2 



2 



?(f) 



^0 I 2 



0.260542 



2^3 77 



„ii3 



= -0.136769 



(59) 
(60) 

(61) 
(62) 



and where we use expressions (||) for xi{t), X2{t) and x^it). 

To calculate the coefficient of restitution one has to use the concept of inverse collision, as was introduced in |3|] 
and discussed in previous chapters of the present study. One obtains the solution of this inverse collision by replacing 
g ^ eg for the initial velocity and a — *■ —a for the dissipative coefficient. In particular, this applies to the maximal 
compression of the inverse collision x™^x — 2;max(5 ^g,ct ^ —a). For consistency one has to require the maximum 
compressions for direct and inverse collision to be equal, i.e. 



(63) 



or using (jSq), 



= g"- [ Vo + ag'-yi + a^g^y2 + a^g^'Vs + 



(64) 



Eq. (^) is in fact an algebraic equation for e^/*^, which may not be generally solved. For this reason we write e as 
an expansion of ag^^^, which is the only combination in which both parameters appear 
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and substitute (^3) into Collecting orders we find 



ag 



1/5 



6 



-^^3 + 25^1 <^2 - Ti^Q ) yo + yiC2 ~ -^2^1 + 2y: 



125 



(65) 



(66) 



The last Eq. (|6^) yields the final result for the coeSi- 
cients: 



C3 



5 yi 

2^2 

15 

4 U2 
95 /yiy 

315 /yi 
32 U2 



-1.153449 



:Cl = 0.798267 



(67) 



15 y2yi ^ 5^3 

4 2/^ 2ya 
105^22/? , 35 ^3^1 
4 



= -0.483582 



+ 



0.285279 



Using {g*) = Cia, we obtain for coefficients ak in 
expansion (|3^): 



(68) 
(69) 
(70) 
(71) 



ai = 1 

02 = C^2/C2 = 3/5 

03 = C^/Cf = 0.315119 

04 = C^/Cf = 0.161167. 



Note that although the general method given in Ap- 
pendix allows to evaluate up to a desired precision all, in 
principle, coefficients C^, it does not provide the closed 
form expression for Ci as the simple approximate ap- 
proach given in the main text does. 
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